% Calculate logit probability for chosen alternatives for each person
%    with multiple choice situations for each person and multiple people,
%    using globals for all inputs except coefficients.
% Original version written by Kenneth Train--first version on July 14, 2006..lastest edits on Sept 24, 2006.
%
% Edits to include discounting by Heather Snell...last edits on Dec 18, 2020

% Simulated Mixed Logit Probability and gradient
% Logit probability is Prod_t { exp(V_*t)/SUM_j[exp(V_jt)] } 
%             = Prod_t { 1 / (1+ Sum_j~=* [exp(V_jt-V-*t)] }
%    where * denotes chosen alternative, j is alternative, t is choice situation.
% Using differences from chosen alternative reduces computation time (with
%    one less alternative), eliminates need to retain and use the dependent
%    variable, and avoids numerical problems when exp(V) is very large or
%    small, since prob is 1/(1+k) which can be evaluated for any k, including
%    infinite k and infinitesimally small k. In contrast. e(V)/sum e(V) can
%    result in numerical "divide by zero" if denominator is sufficiently small
%    and NaN if numerator and denominator are both numerically zero.
% 
% Input f contains the fixed coefficients, and has dimension NFX1.
% Input c contains the random coefficients for each person, and has dimensions NV x NP.
% Either input can be an empty matrix. 
% Output p contains the logit probabilities, which is a row vector of dimension 1xNP.
% Output g contains the gradients of log(p), which is a matrix (NF+NV+NV) x NP
% Code assumes that all GLOBALS already exist.


function [p]=llgrad2nog(f,b,w); 

global NV NF NP NDRAWS NALTMAX NCSMAX NMEM NTAKES XX ANA
global X S XF DR XMAT T
global MDR discountfun FirstRun 

p=zeros(1,NP);

c=trans(b,w,DR);   %Transforms draws into random coefficients c is NV x NP x NMEM
c=permute(c,[2,1,3]); %now NP by NV by NMEM

% This is where discounting occurs and is the most significant difference
% from Train's original codes
dodiscountingnog; %this computes v=exp(utilities not chosen- chosen) and gradient grU 

pp=1./(1+sum(v,1)); %pp is 1 x NCSMAX x NP x NMEM--prob of each choice set
%Back to prob
pp=squeeze(pp); %,NCSMAX,NP,NMEM);     %pp is now NCSMAX x NP x NMEM
pp=prod(pp,1);      %pp is 1xNPxNMEM

%sum over draws
pp=(sum(squeeze(pp),2));       %probs is NPx1
p=pp ./ NDRAWS;
p=p';
p(1,isnan(p))=1;
these=find(p(1,:)==0); 
% have to change zeros to something extremely close to zero
%in order to take the log
p(1,these) = 1e-312;